CALIBRATION OF GYROS WITH TEMPERATURE DEPENDENT SCALE 
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ABSTRACT 

The general problem of gyro calibration can be stated as the estimation of the scale 
factors, misalignments and drift-rate biases of the gyro using the on-orbit sensor 
measurements. These gyro parameters have been traditionally treated as temperature- 
independent in the operational flight dynamics ground systems at NASA Goddard Space 
Flight Center (GSFC), a scenario which has been successfully applied in the gyro 
calibration of a large number of missions. A significant departure from this is the 
Microwave Anisotropy Probe (MAP) mission where, due to the high thermal variations 
expected during the mission phase, it is necessary to model the scale factors as functions 
of temperature. 

This paper addresses the issue of gyro calibration for the MAP gyro model using a 
manufacturer-supplied model of the variation of scale factors with temperature. The 
problem is formulated as a least squares problem and solved using the Levenberg- 
Marquardt algorithm in the MATLAB® library function ’NLSQ’. The algorithm was 
tested on simulated data with Gaussian noise for the quaternions as well as the gyro rates 
and was found to consistently converge close to the true values. Significant improvement 
in accuracy was noticed due to the estimation of the temperature-dependent scale factors 
as against constant scale factors. 

INTRODUCTION 

The general problem of gyro calibration can be stated as the estimation of the scale 
factors, misalignments and the biases of the gyro using the on-orbit sensor measurements 
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(ref.l). Usually, the scale factors of the gyros are deemed to be constants which will be 
estimated a priori and only small, possibly slowly varying corrections to these 
manufacturer-supplied values will be estimated during the mission (ref. 2-5). But in some 
missions such as the Microwave Anisotropy Probe (MAP) (ref. 6), due to the high 
thermal variations expected during the mission, the scale factors are themselves functions 
of temperature. In most of these cases, the functional form of the variation of the gyro 
scale factors will be known a priori from the manufacturer and certain associated 
parameters for their actual values on-orbit are to be estimated as part of the gyro 
calibration exercise. 

The existing gyro calibration utilities used at the NASA Goddard Space Flight Center 
(GSFC) assume that the scale factors are constants. Therefore the objective of the work 
reported in this document is to formulate this temperature-dependent gyro calibration and 
provide a methodology and a tool to solve it. 

PROBLEM FORMULATION 

The problem of temperature-dependent gyro calibration is formulated in this section. 

With a view to applying this methodology to the MAP gyro calibration scenario, the gyro 
model conforming to MAP (ref.7) was used here. The 3x1 vector of gyro rates is given 
by: 

~Sj o oTNj" 

oj( t) — M 0 S2 0 N2 ~ b Q) 

0 0 SjjWjJ 

where 

Si -a lo + auv + a l 2 v 2 + a i 3 v 3 , i=l,2,3 (2) 

are the temperature dependent scale factors, v is the voltage, Ni are the gyro telemetry 
counts, b is the 3x1 bias vector to be estimated, and M is the 3x3 misalignment matrix. 
The scale factors are actually functions of temperature but are provided by the 
manufacturer in terms of gyro thermistor voltage variations. 

The gyro calibration problem now reduces to estimating the 28x1 state vector X 

X=[{a ij },{M u }, b T ,qJf (3) 

where {aij} are the 12 coefficients of the scale factors defined in Equation (2), {Mu} are 
the 9 elements of M, and qo is the 4x1 epoch inertial-to-body quaternion . It is convenient 
to define the 3x4 matrix A whose elements are {a t j}. 

It is possible to reduce the number of parameters in the state vector by 3 by redefining the 
alignment matrix to incorporate the linear part of the scale factor corrections by 
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normalization. We are considering applying this reduction in a future version of the 
algorithm. 

The kinematic equation relating the attitude to the gyro rates is (ref. 8) 
q = jQq (4) 

where q is the attitude quaternion and Q is given by 
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Here co = (ct) x ,a) y ,co z ) is the spacecraft angular velocity measured by the gyros with 
components along the body axes. 

X is estimated by minimizing the error between the attitude quaternions computed from 
sensor measurements and the attitude obtained by integrating the kinematic equation 
above. The closed form solution to the equation is given by (ref. 8): 


q(ti+i) = 
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where 

co = \o 2 ( t t ) + co y 2 ( t t ) + co 2 ( ti )] /2 
At = t. . , — t- 

l+l l 
= Q( a)( t t )) 

q(ti) is the attitude quaternion at time U, and q(t i+ i) is the propagated attitude quaternion at 
time U+i . This closed form solution is used in the present formulation to get the 
propagated attitude quaternions q p (t) for evaluating the objective function matrix for 
minimization using the Levenberg-Marquardt algorithm. The quaternion residuals f(X) 
are now defined by 


f(X) = q P (ti) - q,(tt) 


(7) 


where q t ( t) is the true quaternion at time t obtained using the attitude sensor 
measurements. The problem will then be to find the state vector that minimizes the cost 
function 

L(X) = Zi(fi(X)) 2 (8) 
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with ( f(X ) ) 2 defined as the sum of the squares of the individual components of f(X). 

As will be seen below, the above cost function formulation works well. Note that the 
residuals defined in Equation (7) are different from the traditional attitude residuals 
represented by three small rotations about the true attitude (see for example, ref. 3). This 
traditional approach to the current problem is being investigated. Also, by suitably 
introducing a weight matrix in the cost function, non-uniform weighting can be allowed. 
This is not done in the current formulation and will be attempted in the future version. 

SOLUTION METHODOLOGY 

The Levenberg-Marquardt Method 

The general nonlinear least squares problem can be stated as: 

Find X such that L(Y), defined in equation (8), is a minimum. 

The Levenberg-Marquardt (L-M) (ref. 9) method uses a trust-region approach to shrink 
the step sizes to minimize the cost function at each iteration and the state update is given 
by 


Xnew - X old +ln\ (J T J + D).J t F 
where D is the diagonal matrix given by 
D = A. Diag (J T J ), 

A is a multiplication factor, J is the matrix of first partial derivatives dF/dX, and F is the 
residual vector . 

The general procedure for this method is to set A to 1.0 for the first iteration. If the first 
attempt at an iteration reduces the cost function then A is reduced for the next iteration by 
a factor of 10. If the first step increases the cost function, then A is increased by a factor 
of 5 until the cost function is reduced. The final value of A for the iteration is used for the 
next iteration. Rather than compute the sum of squares of fi(X), ‘NLSQ’ requires the user- 
defined function to compute the matrix-valued function 

F(X) = [f x (X)f 2 (X)f 3 (X) ...f m (X)] T 

Here each fi(X) could be a vector in which case the objective function will be matrix- 
valued. 
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The MATLAB Function ’NLSQ’ 

The MATLAB function ‘NLSQ’ is a general function to solve the non-linear least 
squares problem defined in the previous section. The default algorithm used in this 
function is that of L-M algorithm. The calling parameters for this function are 

[X, OPTIONS, F,J] = NLSQ( ’FUN’, XO, OPTIONS, ’GRADFUN’,P1,P2,..) 

This function starts at the state vector Xo and finds a minimum to the sum of squares of 
the functions described in FUN. FUN is usually an M-file which returns a vector of 
objective functions: F=FUN(X). FUN should return F(X) and not the sum-of-squares 
since the cost function is computed implicitly in the algorithm. 

OPTIONS is a vector of optional parameters to be defined. GRADFUN is an optional 
function which returns the partial derivatives of the functions at a given X. If GRADFUN 
is not supplied, numerical derivatives are computed within the function. 

APPLICATION TO SIMULATED DATA 

General Simulation Procedure 

• Assign true values for the misalignment matrix, M, the scale factor 
coefficients matrix. A, the bias vector, b, and the epoch quaternion qo- 

• Assign start and end times U and tf. 

• Assign m, the number of time tags. 

• Assign a temperature/voltage profile for the simulation time interval [ti,tj\. 

• Assume an angular velocity profile for the spacecraft in [U,tf\. 

• Calculate the scale factor vector S =[A ] V , where A is the 3x4 matrix of scale 
factor coefficients, V=[ 1 v v 2 v 3 ] T , v being the voltage as a known function 
of time t. 

• Calculate the telemetry counts after adding the bias using [S']" 1 M" 1 (N(t)+b) 
where [5] is the diagonal matrix of the scale factor components of S. 

• Starting from qo , get the quaternion history using the recursive relation given 
in Equation (6) above, normalize it and save it as "true" attitudes # t (t). 

Specific Simulation Scenario 

The simulated data spanned 1 hour 15 minutes and consisted of inertial pointing periods 
interspersed with three slews of rate 0.1 deg/sec rotating +/-45 deg about the x, y axes 
and +/-90 deg about z axis. The slews lasted 1 hour totally with a gap of 7 minutes 
between slews. The temperature variation in the gyros in the early launch phase was 
taken to be between 5-10 degrees which corresponds to about 1 to 2 volts variation in 
voltage as per discussion with the MAP Attitude Control System (ACS) engineers. Based 
on this, a sinusoidal voltage variation was assigned with an amplitude of 2 volts and 
period of 1 hour, i. e., 
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v = 2sin 2nt/P 


where t is the time and P, the period, is 1 hour. 

The following set of true gyro parameters and epoch quaternion, denoted by the subscript 
"t", were used: 



0.4945513578201 

0.4945513578201 

0.4945513578201 


.00036598584787 

.00036598584787 

.00036598584787 


.00001662910926 

.00001662910926 

.00001662910926 


-.00000087266463 

-.00000087266463 
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Ixi (T 5 


b, = [0.5, 0.5, 0.5] t deg/hr 
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-1 
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0 

0 
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q 0 = [0.5 0.5 0.5 0.5] t 


The quaternion measurements were corrupted with noise equivalent to 40 arc seconds 
root-sum-squared (RSS) in Euler angles to correspond to the MAP star tracker accuracy 
and the gyro measurements were corrupted with Gaussian noise of 0.016 deg/hr (l-o). 


Estimation Steps 

The core of the estimation procedure is the evaluation of the matrix of residuals F(X ): 

• Offset M,A ,b and qo from the true values to get the initial state vector X 0 

• Calculate the gyro rates using S=A V, co= M S N - b where A is the vector of 
telemetry counts 

• Starting from qo, use Equation (6) above to get the propagated attitude 
quaternions, q p (t), at the sampling intervals and normalize it at every time 
instant. 

• Compute the residual vectors via Equation (7) at times tj, .... t m 

• Update the 28x1 state vector X using the L-M algorithm via the MATLAB 
function ’NLSQ’. 


This procedure is repeated till the convergence criterion of the L-M method is achieved. 
In order to scale the independent variables and make their range of variation uniform 
across the state vector, the components of the bias vector and the coefficients of the scale 
factors are multiplied by 10' 6 before using them in the cost function evaluation. It was 
found that without this, the problem becomes ill-conditioned and the procedure 
terminates. 
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DISCUSSION OF RESULTS 


Numerous estimation runs were made for widely ranging values of the a priori state 
vector Xo. Convergence close to the true value was obtained consistently. The following 
is a typical result: 


b = [ 0.51543 0.5052 0.4954 ]deg/hr 


M = 


-1.00000000000000 0.00000000131127 -0.00000000043832 
0.00000000160429 -1.00000000000000 -0.00000000916929 
0.00000000030366 -0.00000001133398 -1.00000000000000 


A = 


0.44509621842136 

0.44509603800549 

0.44509511127285 


0.00032942235090 

0.00032951201204 

0.00033038158483 


0.00001491460493 - 0.00000080505695 
0.00001499181614 - 0.00000080575904 
0.00001508609330 - 0.00000097751931 



q 0 = [ 0.500001 0.500008 0.49994 0.49995] T 


The biases agree to an accuracy of .015 deg/hr and the products of misalignment and 
scale factor matrices match to an accuracy of 10' 10 , which corresponds to 0.017% error. 
Figure 1 gives the plot of the errors between estimated and true Euler angles and Figure 2 
gives the deviations plot of the components of the products of the M and S matrices 
defined in equation (1) as a function of time (i.e., as temperature varies). 


It was found that the solutions were independent of the initial guesses. Some of the initial 
guesses tried were: 0.9 for all diagonal elements of M, 0.1 for all diagonal elements of M, 
zero for all components of b, 2 deg/hr for all components of b, zero for all scale factor 
coefficients and 1.0 for all scale factor coefficients and different combinations thereof. 
Results differing only in the 12 th or 13 th significant figure were obtained with these initial 
guesses, which is essentially due to the limitation of the machine precision. Runs were 
also made not including measurements during maneuver times and the resulting 
accuracies did not vary significantly. 


Further, estimation was carried out with only constant scale factors in the state vector and 
setting all other coefficients to zero as in the conventional gyro calibration methods. As 
seen in Figure 3, which gives the errors in the Euler angles for this case, as compared to 
those of the temperature-dependent gyro calibration case above, the improvement in 
accuracy is significant due to the enhanced state vector. Note that the asterisks in this 
figure represent the same data as in Figure 1, but on a different scale. 


169 









Fig. 2. Deviations Plot of the Components of the Products of Misalignment and 

Scale Factors 
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Fig 3. Euler Angle Errors Improvement Due to Temperature-Dependent Gyro 

Calibration 


CONCLUSIONS 

A gyro calibration problem with temperature dependent scale factors was formulated. 

The Levenberg-Marquardt algorithm in the MATLAB function NLSQ was used to solve 
the resulting least squares problem. Encouraging results using simulated data, conforming 
to the MAP mission, show the feasibility of applying it for recovering the spacecraft gyro 
scale factors. The advantage of this approach is that no partial derivatives of the cost 
function with respect to the state vector are needed. This helps to augment or remove 
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components from the state vector very easily. Also, since the entire state vector is 
estimated using a single span of data, no operator intervention is needed. 

SYMBOLS 

0 ) Vector of gyro rates 

M Misalignment matrix 

S Temperature-dependent scale factors 

N Gyro telemetry counts 

b Bias vector 

X State vector 

ay Scale factor coefficients 

qo Epoch quaternion 

Q Matrix of gyro rate components 

f(X) Vector of residuals 

q t True quaternion 

< 7 P Propagated quaternion 

Sq Quaternion error 

zlt Time interval 
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